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A method based upon the statistical evaluation of the dif- 
ferentiability of the measure along the trajectory is used to 
identify determinism in high dimensional systems. The results 
show that the method is suitable for discriminating stochastic 
from deterministic systems even if the dimension of the latter 
is as high as 13. The method is shown to succeed in identify- 
ing determinism in electro-encephalogram signals simulated 
by means of a high dimensional system. 



I. INTRODUCTION 

Although numerous methods have been developed in 
recent years addressed to detect determinism in time se- 
ries most of them give wrong or ambiguous re- 
sults when applied to high dimensional systems (wc refer 
to systems with attractor's dimension greater than, say, 
five). However, high dimensional systems are ubiquitous 
in nature, as for example in the case of spatially extended 
systems. Therefore it remains the task to develop a ro- 
bust method capable of detecting deterministic behavior 
in systems with many active degrees of freedom ||. 

We hereby claim that, contrary to some widely ac- 
cepted premises in the chaos community it is pos- 
sible to discriminate determinism from stochasticity in 
systems with a high correlation (or information) dimen- 
sion. Our method is based upon a fundamental property 
of deterministic systems, namely, the differentiability of 
the measure along the trajectory. Our numerical imple- 
mentation of this property is robust enough to uncover 
this topological property in short data sets. 

As a very first step, we should clarify what we mean by 
determinism. Attaching to definitions, saying that a time 
series stems from a deterministic system only if it can 
be produced by a discrete- or continuous-time dynam- 
ics without any stochastic source, would isolate a very 
small class of phenomena. In fact, even artificial time se- 
ries from non-stochastic dynamical systems are typically 
affected by numerical noise, which at best can be con- 
sidered to be in its zero- limit. Unfortunately, numerical 



errors tend to be magnified with the reconstruction pro- 
cess and can mimic the effect of a stochastic feedback. It 
is pertinent to note that Casdagli et al. jl(| have shown 
that univariate time series derived from sufficiently high 
dimensional systems cannot in practice be regarded as 
noise-free. 

On the other hand, the importance of detecting a possi- 
ble deterministic dynamics is not merely academic and it 
is often motivated by the intention of controlling the non- 
stochastic part of the system under investigation. Thus, 
operationally, we prefer to imagine an underlying dynam- 
ics composed by a stochastic part (that is, dynamical 
noise) and a non-stochastic one. In this sense, it's not 
always obvious how to measure the weight of one part 
with respect to the other. However, it seems very nat- 
ural to classify as purely stochastic a system which has 
not residual (non-trivial) dynamics if its stochastic part 
is switched off. For example, this is the case of Eqs. (||) 
and (^) that will be discussed later on. More generally, 
we speak of determinism whenever the dynamics has a 
relevant non-stochastic backbone, being conscious that 
we can decide if this is really non-negligible only when 
we adopt a suitable indicator. In this paper we make 
use of the statistical differentiability of the measure as 
introduced in refs. [ fil| , [L2| . There it is argued that this 
quantity is sensible to the amount of dynamical noise 
intrinsic to the system. Here we explore the two oppo- 
site sides of the scenario depicted above, showing that 
it sharply discriminates purely stochastic systems from 
deterministic dynamics, even when the latter takes place 
on rather high dimensional attractors. 

II. METHODS 

A. Tracking the measure along the trajectory 

To illustrate the point, let's consider a dissipative 
dynamical system described by n-first-order differential 
equations x = F(x). The corresponding flow /* maps a 
"typical" initial condition Xo into x(i) = /*(xo) at time 
t. Once transients are over, the motion settles over the 
attractor A. Defined over this set is the natural measure, 
which, from an operational point of view, can be consid- 
ered as the limiting distribution of almost all starting 
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initial conditions, that is, 



fi(B r (x)) = lim - f l Sr (/ T (yo))dr 



(1) 



for almost all yo in the basin of attraction. In Eq. (|l|) 
£? r (x) indicates the hyper-sphere of radius r centered at 
x and ±B r the associated characteristic function. In prin- 
ciple, with infinitely many data points and aribrarily fine 
sampling, one would consider the limit r — > 0. In practice 
a suitable density estimator is needed. Here we adopt a 
fixed- volume Epanechnikov kernel [p|Jl^ 1, with a com- 
pact support B r and a given finite radius r (see below). 
As shown previously Jn],|l2| , smoothness of this measure 
along the trajectory is a good candidate to quantify de- 
terminism. Numerically, the differentiability of the mea- 
sure along the trajectory was evaluated by means of the 
topological statistics developed by Pecora et al. Q . Ac- 
tually we look at the continuity of the logarithmic deriva- 
tive of the measure (see 12 ). 

We want to emphasize that in looking at this property 
of the measure we avoid one of the most common prob- 
lems found in this kind of algorithms, namely, the lack 
of an acceptable scaling region. In fact, this is the main 
problem encountered when dealing with quantities such 
as the maximal Lyapunov exponents and/or Kolmogorov 
entropy |8) , particularly in the case of high dimensional 
systems (see ref. @). Especially troublesome is the cal- 
culation of the information (Z?i) or correlation (D2) di- 
mensions, in which the search for a reliable scaling region 
is the source of most of the misuse of the algorithm. Con- 
trarily, our algorithm needs a realible estimation of the 
natural measure on the attractor not in a range of scales 
but at a fixed one. There is a lower value of resolution 
given by the minimun average interpoint distance, which 
can be estimated as r m i n ~ N~ 2 /° 2 , and an upper limit 
given by the attractor extent (which we always normal- 
ize to the unit hyper-square) p4[ . Below the lower limit 
there would be statistical fluctuations due to the lack 
of neighboring points, while close to the upper limit we 
must confront edge effects [ p~4|JT5| ] . The latter are particu- 
larly important for high embedding dimensions because, 
as the dimension increases, more points stay near the 
attractor boundary. These bounds are by no means ex- 
act, for instance the expression for the lower limit was 
derived by assuming the data points to be uniformly dis- 
tributed over the attractor, which is not always true even 
in cases as simple as nonlinear oscillators. In practice one 
should look in each particular case for the most appro- 
priate scale to fix. As far as the systems discussed here 
are concerned, we have found that setting r at 10% of 
the attractors' linear extent is a good choice. 



B. Statistical tests of continuity 



nates the signal, more discontinuous it becomes. This is 
the case for example of additive noise, e.g. noise added to 
the signal However, this is by no means a general 

rule. For instance, dynamical noise, that is, noise added 
in the equation of motion, is not expected to affect the 
smoothness of the signal. We overcome this drawback 
by using the distribution of points on the trajectory (or 
the natural measure) as a way to evaluate the degree 
of noise in the system (see 12 1). In order to test the 



mathematical properties of the measure, that is, con- 
tinuity, differentiability, inverse differentiability and in- 
jectivity, we borrow the statistical approach developed 
by Pecora et al. [ fl3| |. Basically, the method is intended 
to evaluate, in terms of probability or confidence levels, 
whether two data sets are related by a mapping hav- 
ing the continuity property: A function / is said to be 
continuous at a point xo if Ve > 0, 36 > such that 
I x - x ||< S =>\\ /(x)) - /(x ) ||< e. The results are 
tested against the null-hypothesis, specifically, the case 
in which no functional relation between points along the 
trajectory and the measure exists. Thus, as done in [ ji3| , 
we calculate 



e C o(e) 
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and 



e c o(e,i) = 1- 



Pj 



(2) 



(3) 



where Pj is the probability that all of the points in the 
5-set, around a given point Xj, of the reconstructed tra- 
jectory, fall in the e-set around 



din ^(xj 
di 



The likelihood 



that this will happen must be relative to the most likely 
event under the null hypothesis, P max (see refs. |]l3| and 
fUH). When @c°{ e i3) ~ 1 we can confidently reject the 
null hypothesis, and assume that there exists a continu- 
ous function. As in the work of Pecora et al. jl3| the e 
scale is relative to the standard deviation of the density 
time series, and thus, e £ [0, 1]. Plots of Q c o(e) versus e 
can be used to quantify the degree of statistical continuity 
of a given function. The typical outcome is a sigmoidal 
curve whose width and slope are affected by the level and 
the type of noise contained in the series 

HUH]- In 

or- 



der to characterize the continuity statistics by means of 
a single parameter we can also calculate, 



/ e c o(e)de 
Jo 



(4) 



The limiting values of 8, namely, and 1, correspond to 
a strongly discontinuous and a fully continuous function, 
respectively. Hereafter we shall refer to 9 as CS (Conti- 
nuity Statistics). 



A naive test to quantify noise in signals is to check 
how smooth they are. As long as more noise contami- 



III. RESULTS 
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A. Generalized Mackey-Glass system 

In order to investigate how our algorithm work on 
high dimensional systems, we use a generalization of the 
Mackey-Glass (MG) equation |17|, a delayed feedback 
system, 



x(t) 



ax(t — to) 

l + X W (t-T ) 



y(t) 



y(t) = -u 2 x(t) - py(t) 



(5) 



where a — 3, p — 1.5, u> = 1 and To = 10. As it is stated 
in ref. |Q the Kaplan- Yorke dimension of this system 
is Dky ~ 13.5, which, according to the Kaplan- Yorke 
conjecture, Dky = D%, for a typical attractor. We have 
made a standard reconstruction analysis over time series 
of up to 16384 data points. Although the optimal time 
delay r should be given by the first minimum of the mu- 
tual information, 210 in sampling units, we have used 
a somewhat smaller value, typically 120. Using larger 
delays in high embedding dimensions would reduce dras- 
tically the number of reconstructed points. The choice of 
this value, however, is supported by the fact that in high 
dimensions, the optimal r seems to be smaller than that 
given by the mutual information criterion , at least in 
the case of chaotic continuous-time systems, as our case. 
Surrogate time series Jl^,^] with the same number of 
points have also been generated, using the routines in the 
TISEAN ftf§ package ("endtoend" and "surrogates"). 

We have also calculated the correlation integral of this 
system. As it is well known, this is the most basic proce- 
dure to discriminate between deterministic and stochas- 
tic behavior. The numerical results for the correlation 
integral and its derivative (that gives the correlation di- 
mension) are shown in Figure 1. The results show that, 
for this high dimensional system, no region with constant 
slope is found. This clearly illustrates the well known 
difficulties inherent to the estimation of the correlation 
dimension in high dimensional systems. 

In Figure 2 we report the values of CS for the general- 
ized MG system (reconstruction from the x-time series) 
when the embedding dimension is varied in the range 2- 
30. In order to increase the robustness of the results, we 
have averaged the CS for six differents time series. The 
high values corresponding to the original series, together 
with the lower values of their surrogates, indicate that 
our method is able to identify determinism in a high di- 
mensional system. We also note that whereas the CS for 
the original series remains almost constant for embed- 
ding dimensions larger than the correlation dimension, 
that for the surrogate series decreases steadily. Finally, 
it is here pertinent to observe that the results of refjlj] 
are referred to the standard MG equation (see Eq. ((L3|) in 
the Appendix) , whose correlation dimension is estimated 
around seven while here we have choosed the gener- 
alization (H) just to check that the methodology works 
with more severe test. However, despite the smaller cor- 
relation dimension, in Fig. 8 of [Jl2fl the CS for the stan- 



dard MG system is lower than the CS for the generalized 
version reported here. In the Appendix we discuss the 
solution to this apparent paradox, through a digression 
on the role of the sampling time on the observed levels 
of CS. 



B. Simple stochastic systems 

As already remarked above, the results just discussed 
have to be contrasted with those obtained for stochastic 
systems. Here we show results for two time series derived 
from a linear (|^) and a nonlinear (^) stochastic processes, 
namely, 



x = 9x(t) + T)(t) 



y=(a- 0.5)13 - y(t) + y/2py{t)ri(t) 



(6) 



(7) 



where rj{t) is a Gaussian noise with standard deviation 
0.1, a = (3 = 1, and 9 = —0.9. Both processes are 
examples of purely stochastic systems, as no dynamics 
is left in the noise free limit, and exhibit (l/f a ) power 
law spectra, a finite correlation dimension (« 2) and a 
converging Kolmogorov entropy p2[ . Note that this is 
a typical case in which the correlation dimension fails in 
identifying stochasticity. 

Attractor reconstruction was carried out on time se- 
ries with up to 8192 data points and a time delay of 100 
(somewhat less than the first minimum of the mutual in- 
formation which in this case lies approximately at 130). 
As it is clearly seen in Figure 3, there is no difference 
between the CS for the orginal time series and its surro- 
gate, both showing very low values which indicate a low 
differentiability (a signature of stochasticity as discussed 
in p|). 



C. Model of Electro-encephalogram signals 

A paradigmatic example of high dimensional systems 
is found in physiology, namely, Electro-encephalogram 
(EEG) signals. These are in fact the result of a sum over 
a large number of neuronal potentials. There are many 
studies that claim deterministic behavior in EEG dynam- 
ics, most of them based on the calculation of the corre- 
lation dimension. However, recent analyses have pointed 
out many technical problems related to those studies (see 
Q and references therein) which throw serious doubts on 
the above conclusion. We have applied our method to the 
analysis of a EEG-like signal |)| generated by the follow- 
ing set of nonlinear coupled equations, 



x\ = x 2 



x 5 — 25 . 
x-i = smuit + 3xjsmuJ2t 

+XnSinijJ2^ — 3\X6\X2 — XgXi 



(8a) 



(8b) 
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with uji—30, ^2=65 and W3=80; 

± 3 = <j(x4 - x 3 ) (9a) 

±4 = — .T3.T5 + r.T3 — X4 (9b) 

±5 = X3X4 — bx$ (9c) 
where cr=10, r=28 and 6=8/3 (Lorenz system). 

±6 = £7 (10a) 

£7 = — fcx7 — £g + Boost (10b) 
where fc=0.1 and B=12 (Ueda equations). 

± 8 = x 9 (11a) 

±9 = ~Sxg + 2^8(1 — £g) + fcOSLdt (lib) 

where <5=0. 15, F=0.15 and w=0.8 (two-well potential 
Dufhng-Holmes attractor) . 

iio = -(#11 + £12) (12a) 
in = £10 + «n (12b) 

X12 ^ a + x 12 {x w - fi) (12c) 

where a=0.15 and /i=10 (Rossler attractor). This intri- 
cate system has a correlation dimension around 9 and has 
failed to pass the method originally proposed by Salvino 
and Cawley 

We have integrated the whole system of || using a 
fourth-order Runge-Kutta algorithm with a fixed step of 
0.001. The main two coordinates, x\ and X2, behave in 
a markedly different way. Coordinate x 2 shows a rather 
standard behavior with the first zero of the autocorrela- 
tion function lying at ~ 27 (in units of sampling time). 
Instead, the x\ time series exhibits a power spectrum 
that seems to be of the l// 2 type and a first zero of 
the autocorrelation that strongly depends on the series 
length N. We found 1.3 x 10 4 for N = 10 5 and 9.9 x 10 4 
for N — 10 6 (both in units of a larger time step, namely 
0.01, that was used to analyze longer time series). These 
results may reveal a non-stationary behavior of x\. In 
fact this coordinate shows a saw-tooth shape with a very 
long wavelength and, as a consequence, it may appear 
to increase linearly over rather long time intervals. As 
the analysis of || was done on the X\ coordinate it may 
actually be the reason why those authors failed in de- 
tecting determinism in this system. Note that the lack 
of stationarity (a property related to ergodicity) would 
imply the failure of our method. A proper estimation 



of the (reconstructed) measure, requires that space aver- 
ages correspond to time averages, which in fact would be 
impossible if the time series is not stationary. This is the 
reason why we cannot estimate confidently the measure 
along the trajectory using x\ in the above system. How- 
ever, as x\ = X2, this allow us a better estimate from 
the X2 coordinate (since upon differentiating jl|] the saw- 
tooth shape of the xi-time series coordinate results in 
very small offsets). 

Figure 4 shows the CS for both the X2 time series and 
its surrogate. We have averaged over 8 differents real- 
izations and their corresponding surrogates. The results 
correspond to time series having 4096 data points and a 
time delay of 27 (in this case the first zero of the auto- 
correlation function and the first minimum of the mutual 
information almost coincide). The difference between the 
CS for the original series and that for its surrogate is no- 
ticeable, even on such short time series. This, along with 
the high value of CS obtained for the original series |fL2| , 
reveals an essentially deterministic origin of the simulated 
EEG signal. 

IV. DISCUSSION 

The examples discussed in the previous section indi- 
cate that the analysis of the CS with respect to the em- 
bedding dimension reveals useful informations in order 
to discriminate stochasticity from determinism. Despite 
the high dimensionality of the deterministic systems we 
notice that the distinction is feasible almost from the 
start, that is, from dimensions lower than the attractor 
dimension. This is an interesting (though somehow "for- 
tunate" ) outcome since, for such low embedding dimen- 
sions, transverse self-intersections of the orbit must occur 
(or equivalently the existence of false nearest neighbors). 
However, as the embedding dimension increases, and es- 
pecially beyond the correlation dimension of the system, 
the CS of the original signal seems to stabilize around 
a constant value. Contrarily, the CS for the surrogated 
series steadily decreases for embedding dimensions larger 
than the correlation dimension (compare Figures 2 and 
4). Finally, the CS for the purely stochastic processes 
(Figure 3) attains rather low values and is monotonically 
decreasing. 

The small error bars in the CS, both in Figure 2 and 
4, indicate that the smoothness of the measure in de- 
terministic systems shows almost no dependence on ini- 
tial conditions. The absolute value of the CS for deter- 
ministic signals may slightly depend on several factors, 
namely, number of data points, sampling rate, ball size 
in the measure estimation process, etc. But, for most 
deterministic time series belonging to different trajecto- 
ries, the computed CS seems to be a robust measure of 
the property we want to identify. On the other hand, the 
CS calculated for the reconstructed measure in the case 
of the surrogate time series shows a large variability. In 
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noting this we are not only referring to the large error 
bars of Figures 2 and 4 characteristic of most noisy mag- 
nitudes, but rather to the largely different CS found for 
the systems of those two Figures. In particular we note 
that, while the CS for the original series are similar and 
higher than 0.9, those for the surrogated series are, higher 
than 0.6 in the case of the generalized Mackey- Glass (Fig. 
2) and smaller than 0.2 for the EEG signal (Fig. 4). Sur- 
rogate time series generated by the IAAFT method [^0| , 
as those used here, are in fact, realizations of a Gaussian 
linear stochastic process (possibly passed through a non- 
linear invertible measurement function). By definition, 
they share the same correlation structure of the origi- 
nal time series and the same probability density function 
(histogram). In passing, notice that in refs. 0,^2| the 
results correspond to the simpler AAFT method. Al- 
though correlation structure is in some way related to 
continuity, both are very different properties of the sig- 
nal. A strict application of the surrogate methodology in 
our case would requiere to generate surrogates time series 
with the same continuity characteristics of the original 
signal. Moreover, and going back to the issue of the ab- 
solute value of the CS, we should keep in mind that the 
statistics of Pecora, Carroll and Heagy |l3| are in turn a 
measure of a certain null-hypotesis. Strictly speaking we 
are dealing with the superposition of two null-hypotesis 
and this drawback doesn't allow a clear interpretation of 
the CS values coming from surrogate data. It is likely 
that a different surrogation method is required for a safe 
application of our method to the analysis of surrogate 
series derived from nonlinear systems. We are actually 
investigating this issue in the light of constrained ran- 
domization l2~3]. 



V. CONCLUDING REMARKS 

Summarizing, we have shown numerically that it is 
possible to discriminate determinism from stochasticity 
in high dimensional systems. A local property of the 
system, as the differentiablity of the measure along the 
trajectory, prevents us from regarding as "random" those 
systems that actually are deterministic (a task in which 
other methods have failed) . The fact that for sufficiently 
fine sampling and large embedding dimensions the pa- 
rameter CS for deterministic systems is always above 0.9, 
opens the possibility of applying our methods to the anal- 
ysis of more complicated series such as the experimental 
ones without resorting to comparison with the CS for 
surrogate series. A lot of work still remains to establish 
the method as a truly practical tool. For instance, its 
limits of applicability have to be identified, particularly 
regarding the relation between the minimum number of 
data points in the series and the actual dimension of the 
system. Real applications would require to extract in- 
formation about the existence of many more (hundred, 
thousand, etc.) active degrees of freedom of the under- 



lying system. Fortunately, as the results here discussed 
seem to indicate, this does not necessarily require large 
embedding dimensions. 
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VI. APPENDIX: HOW RESULTS DEPEND ON 
SAMPLING 

One of the very prcliminar steps in time scries analysis 
in to ensure that one is dealing with a good sampling 
time. The precise meaning of this statement depends 
somehow on the specific problem at hand. Nonetheless, 
what is normally done is to check that the sampling time 
is much smaller than the smallest time scale, i m i n , present 
in the system. This could be a (quasi) period, an auto- 
correlation time, an inverse of the maximum Lyapunov 
exponent and so on. However, there are practical limi- 
tations to this rule. Experimentally, the resolution time 
is not always freely adjustable. Numerically, one should 
cope with the limited performances of the machine. In 
both cases, one has to consider also the problem of storing 
large amount of data. This is fundamental when one re- 
alizes that the system possesses intrinsic long time scales, 
which have to be covered in order to analyze stationary 
time series. The compromise between series length and 
stationarity is often achieved by making a re-sampling of 
the series itself. However, this decimation procedure is 
not always innocuous. This Appendix is devoted to study 
how our basic index, the CS, is affected when one takes 
different sampling time. We take the following three ex- 
amples. First, time series built from the ^-coordinate of 
the Lorenz system (x^ in Eqs. (g)). Second, the Mackey- 
Glass delayed differential equation: 



ax(t — To) 

l + X 10 {t-T ) 



bx(t) , 



(13) 



choosing a — 0.2, b = 0.1 and To = 100. Third, the purely 
(nonlinear) stochastic system of Eq. ([?]). In Fig. 5 we 
report the corresponding CS's versus the embedding di- 
mension, for various values of the decimation time. The 
resulting series lengths are 10000 for the Lorenz and MG 
systems and 20000 for the stochastic process. In every 
case, the delay time of the reconstruction was choosen as 
the minimum of the corresponding mutual information. 
As far as the Lorenz and MG systems are concerned, one 
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sees what it is somehow expected for a discretized con- 
tinuous function, that is, the coarser is the sampling the 
lower is the continuity level. Hence, when applying this 
kind of methods, one should keep in mind that to some 
extent the CS may be lowered due to a not-so-good sam- 
pling. In particular from Figs. 2 and 5 one sees that is in 
fact the or igin of the discrepance, anticipated in Subsec- 
between the MG data shown here and those 
Nevertheless, an important feature must be 
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tion 
of ref. 

pointed out, namely that this "runoff" is not arbitrary. 
The CS for the determinstic systems doesn't fall down to 
the values corresponding to the stochastic system which 
in turn appear to be essentially independent of the sam- 
pling time. It is interesting to speculate to what extent 
this last property can be exploited to build up a tech- 
nique based on re-sampling, where the stochastic signals 
are identified just as those which are "squeezed down" to 
rather low values of CS, independently of their resolution. 
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FIG. 1. Log-Log plot of the correlation integral, with em- 
bedding dimensions in the range 10 to 40, for the y-coordinate 
of system (^). 16384 data points has been used. Upper lines 
show the derivative of the correlation integral for embedding 
dimensions of 20, 22, 24, 28, 30, 32, 34 and 36. Horizontal 
dashed line shows the actual information dimension of the 
system. 



1.0 



0.8 



0.6 



CO 
O 



0.4 



0.2 



0.0 



10 20 
Embedding Dimension 



30 



G 



FIG. 2. CS for the system (|5|) (filled circles), averaged over 
six differents time series of 16384 data points each. Error bars 
(standard deviation) are also shown. The {/-coordinate has 
been used. The empty symbols correspond to their respective 
surrogates. 



FIG. 4. CS for the EEG-like signal. The results correspond 
to the X2 coordinate of the system. 4096 data points were used 
in each of the eight time series averaged. The standard devi- 
ation (error bars) is also shown. Empty symbols correspond 
to their respective surrogates. 
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FIG. 3. CS for random signals, a) Solid circles correspond 
to Eq. Q and empty circles are the corresponding surro- 
gates (differences are not visible as the points almost coin- 
cide). Eight differents realizations of 8192 each have been 
used, b) Idem as a) but for Eq. (^). 
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FIG. 5. Dependence of the CS on the sampling time for dif- 
ferent continuous-time systems. Upper panel: Lorenz system; 
Central panel: Standard Mackey-Glass equation (Eq. (|l3|)); 
Lower panel: Random system of Eq. (Q) . The symbol dt indi- 
cates the time step used in the numerical integration scheme 
(fourth-order Runge-Kutta for the first two and Euler's for 
the random system) while n corresponds to the re-sampling 
step. 
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